Mathematical analysis of stochastic 
models for tumor-immune systems 



O. Chi§*, D. Opri§* 



* Euro University "Dragan", Lugoj, Romania 
** Faculty of Mathematics and Informatics, West University of Timi§oara, Romania 
E-mail: chisoana@yahoo.com, opris@math.uvt.ro 

Abstract: In this paper we investigate some stochastic models for tumor- 
immune systems. To describe these models we used a Wiener process, as the 
noise has a stabilization effect. Their dynamics are studied in terms of stochas- 
tic stability in the equilibrium points, by constructing the Lyapunov exponent, 
depending on the parameters that describe the model. We have studied and 
analyzed a Kuznetsov- Taylor like stochastic model and a Bell stochastic model 
for tumor-immune systems. These stochastic models are studied from stability 
point of view and they were represented using the Euler second order scheme. 
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1 Introduction 

Cancer is a disease that may affect people at all ages. It causes about 13% 
of all human deaths. In the prognosis of cancer patients, it should be taken into 
consideration the type of cancer and the stage of the disease. Cancers may be 
treated or cured, depending on the specific type, location, and stage. We may 
say that surgery and chemotherapies play an important role in treating cancer, 
but they do not represent a cure. What it is needed is a successful treatment 
strategies, one of these strategies is investigated through immunotherapy [P3j . 
by defining a model of differential equations that represents the interaction 



1 



2 



between effector cells and tumor cells. This idea of immunotherapy is promis- 
ing, but controversial from the point of view of the results obtained in medical 
investigations. 

Stochastic modelling plays an important role in many branches of science. 
Because in practical situations we confront with instability and perturbations, 
we will express our mathematical models using white noise, represented by 
brownian motion. We will study stochastic dynamical systems that are used 
in medicine, in describing a tumor behavior. Cancer tumor may be destroyed 
using some treatments, but a regression of the disease may appear. So, we need 
not only preventative measures, but also more successful treatment strategies. 
Efforts along these lines are now being investigated through immunotherapy 
(0, [20], [22])- A simulating model is described by the existence of tumor 
free equilibrium. A tumor size may tend to +oo, depending on the parameters 
of the model, and may exist a "small tumor size" equilibrium, which coexists 
with the tumor free equilibrium [TO] . 

This tumor-immune study, from theoretical point of view, has been done for 
two cell populations: effector cells and tumor cells. It was predicted a thresh- 
old above which there is uncontrollable tumor growth, and below which the 
disease is attenuated with periodic exacerbations occurring every 3-4 months. 
There was also shown that the model does have stable spirals, but the Dulac- 
Bendixson criterion demonstrates that there are no stable closed orbits. It is 
consider ODE's for the populations of immune and tumor cells and it is shown 
that survival increases if the immune system is stimulated, but in some cases 
an increase in effector cells increases the chance of tumor survival. 

In the last years, stochastic growth models for cancer cells were developed. 
These models simulate the way tumors evolve with respect to a certain therapy, 
but also they show the interactions between tumor cells and immune cells. We 
mention the papers of W.Y. Tan and C.W. Chen [TH], N. Komarova, G. Albano 
and V.Giorno pQ, L. Ferrante, S. Bompadre, L. Possati and L. Leone [6], A. 
Boondirek Y. Lenbury, J. Wong-Ekkabut, W. Triampo, I.M. Tang, P. Picha 

0- 

Our goal in this paper is to construct stochastic models and to analyze 
their behavior around the equilibrium point. In these points stability is stud- 
ied by analyzing the Lyapunov exponent, depending of the parameters of the 
models. Numerical simulations are done using a deterministic algorithm with 
an ergodic invariant measure. In this paper, the authors studied and analyzed 
two stochastic models. In Section 2, we considered a Kuznetsov and Taylor 
stochastic model. Beginning from the classical one, we have studied the case 
of positive immune response. We gave the stochastic model and we analyzed 
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it in the equilibrium points. Numerical simulations for this new model are 
presented in Section 2.1. In Section 3 we presented a general family of tumor- 
immune stochastic systems and from this general representation, we analyzed 
Bell model. We wrote this model as a stochastic model, using Annexe 1, and 
we discussed its behavior around the equilibrium points. Numerical simula- 
tions were done using the software Maple 12 and we implemented the second 
order Euler scheme for a representation of the discussed stochastic models, 
described in Annexe 2. 



2 Kuznetsov and Taylor stochastic model 

We will begin our study from the model of Kuznetsov and Taylor [13] . This 
model describes the response of effector cells to the growth of tumor cells and 
takes into consideration the penetration of tumor cells by effector cells, that 
causes the interaction of effector cells. This model can be represented in the 
following way: 

J ±(t) = ai - a 2 x(t) + a 3 x(t)y(t), , . 

\ y(t)=b iy (t)(l-b 2 y(t))-x(t)y(t), [L) 

where initial conditions are x(0) = xq > 0, y(0) = yo > and 03 is the immune 
response to the appearance of the tumor cells. 

In this paper we consider the case of 03 > 0, that means that immune re- 
sponse is positive. For the equilibrium states Pi and P 2 , we study the asymp- 
totic behavior with respect to the parameter a\ in (TjQ). For bia 2 < ai, the 
system (TjQ) has the equilibrium states P\{xi,yi) and ^2(^2, 2/2), with 

xi = -,Vi = 0, (2) 
a 2 

x 2 = (61 (a 3 - b 2 a 2 ) + v / A)/(2a 3 ), y 2 = (h(a 3 + b 2 a 2 ) - v / A/(26 1 6 2 a 3 ) (3) 
where A = b\{b 2 a 2 — a 3 ) 2 + 4bib 2 aia 3 . 

In [13] it is shown that there is an a 10 such that if a% < a w , the equilibrium 
state Pi is asymptotical stable, for ai > a w the equilibrium state Pi is unstable 
and if a\ < a\o the equilibrium state P2 is unstable and for a\ > aio the 
equilibrium state P2 is asymptotical stable. 

In the following, we associate a stochastic system of differential equations 
to the classical system of differential equations ([T]). 



Let us consider (fi, £F t > , O 5 ) a filtered probability space and (W(t)) t >o a 
standard Wiener process adapted to the filtration (J) t > . Let {X(t) = (x(t),y(t))} 
be a stochastic process. 

The system of It 6 equations associated to system (0Q) is given by 



x 



(t) = x + / *(ai - a 2 x(s) + a 3 x(s)y(s))ds + gi{x{s),y{s))dW{s) 1 

V(t) = y + j*{{b iy {s){l - b 2 y{s)) - x{s)y{s))ds + g 2 {x{s),y{s))dW{s), 

(4) 

where the first integral is a Riemann integral, and the second one is an Ito 
integral. {W(t)} t > is a Wiener process [To] . 

The functions gi(x(t), y{t)) and g 2 (x(t),y(t)) are given in the case when we 
are working in the equilibrium state. In P\ those functions have the following 
form 

gi(x(t),y(t)) = b u x(t) + b l2 y(t) + c n , 

(5) 

g 2 (x(t),y(t)) = b 21 x(t) + b 22 y(t) + c 2 i, 

where 

cn = -bnxi - b 12 y u c 21 = -b 2l x x - b 22 y x . (6) 

In the equilibrium state P 2 , the functions gi(x(t),y(t)) and g 2 (x(t),y(t)) 
are given by 

gi(x(t),y(t)) = b n x(t) + b 12 y(t) + c 12 , 

(7) 

92(x(t),y(t)) = b 21 x{t) + b 22 y{t) + c 22 , 

where 

C12 = -611X2 - b 12 y 2 , c 22 = -621^2 - b 22 y 2 . (8) 

The functions gi(x(t),y(t)) and g 2 (x(t),y(t)) represent the volatilisations 
of the stochastic equations and they are the therapy test functions. 



2.1 The analysis of SDE ((H). Numerical simulation. 

Using the formulae from Annexe 1, Annexe 2, and Maple 12 software, 
we get the following results, illustrated in the figures below. For numerical 
simulations, we use the following values for the parameters of the system (jlj): 

ai = 0.1181, a 2 = 0.3747, a 3 = 0.01184, h = 1.636, b 2 = 0.002. 
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The matrices A and B are given, in the equilibrium point Pi{^, 0) by 



A 



-a 2 + a 3 yi a 3 x 1 
-yt 61 - 26 2 2/i - x\ 



B 



10 -2 
2 10 



In a similar way, matrices A and B are defined in the other equilibrium 
point 

-hfactz - a 3 ) + V&) (61(6202 + 03)-^)- 



Pi 



2a, 



26162O3 



with A = 6f (6 2 a 2 — as) 2 + 461620103. 



Using the second order Euler scheme for the ODE system ([T]), respectively 
SDE system 01]), we get the following orbits. 
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Fig 1: (n,x(n)) in Pi for ODE © 
optimal behavior of tumor cells 
for ODEQ in Pi 




Fig 3: (n, y{n)) in Pi for ODE Q 
optimal behavior of effector cells for 
for ODEQ in P 1 



0.316- 
0.314- 
0.312- 
0.310- 
0.308- 
0.306- 
0.304- 
0.302- 
0.300- 
0.298- 



Fig 2: (n,x(n,u>)) in Pi for SDE gjl 
optimal behavior of tumor cells 
for ODEgJ in Pi 
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Fig 4: (n,y(n,u)) in Pi for SDE @) 
optimal behavior of effector cells 
for ODElgJ in Pi 
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Fig 5: (x(n),y(n)) in Pi for ODE (T) 
optimal behavior of tumor cells 
vs effector cells for ODEQ in Pi 




Fig 7: {n, x(n)) in P 2 
optimal behavior of tumor cells 
for ODEjTJ) in P 2 



Fig 6: (x(n,w),y(n,w)) in P 1 for SDE © 
optimal behavior of tumor cells 
vs effector cells for ODE© in Pi 
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Fig 8: (n, x(n,w)) in P2 
optimal behavior of tumor cells 
for ODE® in P 2 



Fig 9: (n, y(n)) in P2 
optimal behavior of effector cells 
for ODEjlJ in P 2 



Fig 10: (n,y(n,w)) in P2 
optimal behavior of effector cells 
for ODE© in P 2 
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Fig 11: (x(n),y(n)) in P2 Fig 12: (x(n,u)),y(n,u))) in P2 

optimal behavior of tumor cells optimal behavior of tumor cells 

vs effector cells for ODEQ in P2 vs effector cells for ODEQ in P2 




Fig 13: (a,X(a)) in Pi Fig 14: (a,X(a)) in P 2 



The Lyapunov exponent, for the equilibrium point P\ is negative, so P\ 
is asymptotically stable for each a G K. For the equilibrium point P2, it is 
asymptotically stable for all values of a from the interval (—1.8,1.8), that 
means that P 2 is unstable for all a G (—00, —1.8) U (1.8, 00). 

3 A general family of tumor- immune stochas- 
tic systems 

A Volterra-like model was proposed in [T7], for the interaction between a 
population of tumor cells (whose number is denoted by x) and a population 
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of lymphocyte cells (y), and it is given by 

x(t) = ax(t) - bx(t)y(t) 



y(t) = dx(t)y(t)-fy(t)-kx(t), (9) 

where the tumor cells are supposed to be in exponential growth (which is, how- 
ever, a good approximation only for the initial phases of the growth) and the 
presence of tumor cells implies a decrease of the "input rate" of lymphocytes. 

A general representation for such models can be considered in the form 
given by d'Onofrio in [5]: 

x(t) = fi(x(t),y(t)), y{t) = f 2 (x(t),y(t)), 
x(0) = x , y(0) = yo, 1 ] 

where x is the number of tumor cells, y the number of effector cells of immune 
system and 

fi(x(t),y(t)) = x(t)(hi(x(t)) - h 2 (x(t))y(t)), 

(11) 

f 2 (x(t), y(t)) = (h(x(t)) - h 4 (x(t)))y(t) + h(x(t)). 

The functions hi, h 2 , h 3 , h±, h 5 are given such that the system ffTUl) admits the 
equilibrium points Pi(xi,yi), with x\ — 0, yx > 0, and P 2 {x 2 ,y 2 ), with x 2 7^ 
0, y 2 > 0. 

Particular cases, that will be discussed here, are the following: 

Volterra model |21j ifhi(x(t)) = ai, h 2 (x(t)) = a 2 x(t), hs(x(t)) = b 3 x(t), h^(x(t)) = 
b 2 and h${x{t)) = —bix(t); 

Bell model |3J hi(x(t)) = aix(t), h 2 (x(t)) = a 2 x{t), h s (x(t)) = bix(t), h 4 (x(t)) = 
b 3 and h 5 (x(t)) = —b 2 x(t) + 6 4 ; 

Stepanova model [18J with hi(x(t)) = ai, h 2 (x(t)) = 1, h 3 (x(t)) = b\x(t), hi(x(t)) = 
b and h^(x{t)) = —b 2 x{t) + 6 4 ; 

Vladar-Gonzalez model |20j if in ffTUl) we consider /ti (x(t)) = log(K/x(t)), h 2 (x(t)) 
1, h 3 (x(t)) = hx(t), h 4 (x(t)) = b 2 + b 3 x 2 (t) and h 5 (x(t)) = 1; 

Exponential model [22] if in (JTUJ) we consider hi(x(t)) = 1, h 2 (x(t)) = 1, 
h 3 (x(t)) = bix(t), h 4 (x{t)) = b 2 + b 3 x 2 (t), and h 5 (x(t)) = 1; 

Logistic model |14| if in (JTUJ) we consider hi(x(t)) = 1 — h 2 (x(t)) = 1, 
h 3 (x(t)) = bix(t), h 4 (x(t)) = b 2 + b 3 x 2 (t), and h 5 (x(t)) = 1. 
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For a considered filtered probability space (Q,3 r t >o,y) and a standard 
Wiener process (W(t)) t >o, we consider the stochastic process in two dimen- 
sional space (90t>o- 

The system of Ito equations associated to system ( ITOl) is given, in the 
equilibrium point P(xo,yo), by 

x(t) = x + Ux(s)(hi(x(s)) - h 2 {x(s))y(s)]ds + f* gi(x(s),y(s))dW(s), 
y(t) = Vo + Jo [(h 3 (x(s)) - h 4 (x(s)))y(s) + h 5 (x(s))]ds+ 
+ J*g 2 (x(s),y(s))dW(s), 

(12) 

where the first integral is a Riemann integral, and the second one is an Ito 
integral. {W(t)} t >o is a Wiener process [16]. 

The functions gi(x(t), y(t)) and g2(x(t),y(t)) are given in the case when we 
are working in the equilibrium state P e , and they are given by 

gi{x(t), y(t)) = b n x(t) + b 12 y(t) + c le , 

(13) 

g 2 {x{t), y{t)) = b 21 x(t) + b 22 y(t) + c 2e , 

where 

c ie = -b a x e - b i2 y e , i = l,2, (14) 

and bij e M, i,j = 1, 2. 



3.1 Analysis of Bell model. Numerical simulations. 

Following the algorithm for determining the Lyapunov exponent (Al) and 
the description of the second order Euler scheme (A2) in Maple 12 software, 
we get the following results, illustrated in the figures below. For numerical 
simulations we use the following values of parameters: 

a t = 2.5, a 2 = 1, b x = 1, b 2 = 0.4, b 3 = 0.95, 6 4 = 2. 

The matrices A and B are given, in the equilibrium point P\ by 

A = (~ a2yi + ai ~ a 2 x i \ b = ( a 

' "\-b 2 + b iyi hxi-hj' \P «J' 

with a = a G M, f3 = —2. In a similar way the matrices A and B are defined 
in the equilibrium point P 2 ( a ^~ a ^ 9±] . 
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Fig 15: (n, x(n)) in Pi for ODE JTTJ 
optimal behavior of tumor cells 
for ODEJTT) in Pi 



100 200 300 400 500 600 700 800 

Fig 16: (n,x(n,uj)) in Pi for SDE (IS) 
optimal behavior of tumor cells 
for ODE(T2) in Pj 
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Fig 17: in, y(n)) in Pi for ODE JTTJ 
optimal behavior of effector cells 
for ODEl(TlJl in Pi 
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Fig 18: (n,y(n,uj)) in P x for SDE (12) 
optimal behavior of effector cells 
for ODE(l2) in Pi 
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Fig 19: (x(n),y(n)) in Pi for ODE (TT) 
optimal behavior of tumor cells 
vs effector cells for ODEJTT) in Pi 



Fig 20: (x(n,uj),y(n,uj)) in Pi for SDE (12) 
optimal behavior of tumor cells 
vs effector cells for ODE(l2) in Pi 





Fig 25: (x(n),y(n)) in P 2 for ODE {TTJ 
optimal behavior of tumor cells 
vs effector cells for ODEJTTJ in P 2 



Fig 26: (x(n,uj),y(n,uj)) in P 2 for SDE JT2} 
optimal behavior of tumor cells 
vs effector cells for in P 2 
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The Lyapunov exponent variation, with 6 n = a a variable parameter, 
is given in Figure 27 for the equilibrium point P 1; and in Figure 28 for the 
equilibrium point P2- 




Fig 27: (a, \(a)) in Pi Fig 28: (a, A(a)) in P 2 



From the figures above, the equilibrium points Pi and P2 are asymptotically 
stable for all a such that the Lyapunov exponents A (a) < 0, and unstable 
otherwise. So, Pi is asymptotically stable for a G (—00,— 1.78) U (2.02, 00) 
and P2 is asymptotically stable for a G (—00, —1.62) U (1.88, 00). 

4 Conclusions 

As considered in this paper, we used established conceptual models, but it 
is also very important to consider the model through all its aspects, as we have 
done in this case by imposing the positivity of its solutions. Even if the initial 
model violates the positivity rule, it is valuable because it may be read as a 
model which takes into account a disease-induced depression in the influx of 
lymphocytes. Then, instead of proposing another specific model, we preferred 
to add this new feature to a family of equations, and so, in particular to our 
models chosen for study. 

We have focused on two important tumor-immune systems, presented from 
stochastic point of view: a Kuznetsov- Taylor model and Bell model, that be- 
longs to a general family of tumor-immune stochastic systems. We have deter- 
mined the equilibrium points and we have calculated the Lyapunov exponents. 
A computable algorithm is presented in Al. These exponents help us to decide 
whether the stochastic model is stable or not. For numerical simulations we 



13 



have used the Euler scheme presented in detail in A2 and the implementation 
of this algorithm was done in Maple 12. In a similar way other models given by 
pip can be studied. The model given by the SDE (fl2l) allows the control of the 
model given by ODE (CQ) with a stochastic process. This model is dependent 
on initial conditions. These are very difficult to find for a concrete case, that 
is why it is quite impossible to plan an anticancer therapy based only on this 
method. This is the ony disadvantage for the immunotherapy. 

In our further work, we will consider the tumor-immune model with delay, 
and also another technique used for a successful therapy, using synchronization 
of the coupled tumor-immune model of repressilators in tumor cells aggrega- 
tions. 



Annexe 

Al Lyapunov exponents and stability in stochastic 2- 
dimensional structures. 

The behavior of a deterministic dynamical system which is disturbed by 
noise may be modelled by a stochastic differential equation (SDE). In many 
practical situations, perturbations are generated by wind, rough surfaces or 
turbulent layers are expressed in terms of white noise, modelled by brownian 
motion. The stochastic stability has been introduced by Bertram and Sarachik 
[T2] and is characterized by the negativeness of Lyapunov exponents. But it 
is not possible to determine this exponents explicitly. Many numerical ap- 
proaches have been proposed, which generally used simulations of stochastic 
trajectories. 

Let (Q, J, IP) a probability space. It is assumed that the a— algebra 9^(i > 
0) such that 

^CJ ( C J, Vs < t, s,t E I, 

where I — [0,T], T e (0,oo). 

Let {x(t) = (xi(t), a; 2 (t))}t>o be a stochastic process. The system of Ito 
equations 

dxi(t, u) = fi(x(t, u))dt + gi(x(t, ou))dW(t, u), i = 1, 2, (15) 
with initial condition x(0) = Xq is interpreted in the sense that 

Xi(t,uj) = x i0 (t,uj)+ / fi(x(s,uj))ds+ / gi(x(s,uj))dW(s,uj), i = 1,2, (16) 
Jo Jo 
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for almost all uj G Vt and for each t > 0, where fi(x) is a drift function, 
Qi{x) is a diffusion function, f* fi(x(s))ds, i — 1,2 is a Riemann integral and 
Jo 9i( x ( s ))dW(s), i = 1,2 is an Ito integral. It is assumed that fi and g i7 i = 
1,2 satisfy the conditions of existence of solutions for this SDE with initial 
conditions x(0) = a$ G W 1 . 

Let so = (stab ^02) G M 2 be a solution of the system 

/i(x ) = 0, i = l,2. 
The functions ^ are chosen such that 

9i{x ) = 0, i = 1,2. 
In the following, we will consider 



(17) 



where 6y G 



z = 1,2. 

The liniarized of system ( 1T61) in xq is given by 



X(t) = [ AX(s)ds+ [ BX(s)dW(s) 
Jo Jo 



where 



X(t) 



x(t, uj) 
y(t,u) 



A = 


a n 


a\2 


, B = 


fen 612 










&21 ^22 



dfi 



dxj 



dgi 



8xa 



(18) 



(19) 

(20) 
(21) 

XO UJjj XO 

The Oseledec multiplicative ergodic theorem [15] asserts the existence of 
two non-random Lyapunov exponents A2 < Ai = A. The top Lyapunov expo- 
nent is given by 

A = lim - sup log ^x(t) 2 + y(t) 2 . (22) 

i-+oo t 

Applying the change of coordinates 

x{t) = r{t) cos0(*), x{t) = r{t) sin0(t), 
by writing the Ito formula for 

h i{x, y) = - log(x 2 + y 2 ) = log(r), 



h 2 (x, y) = arctan ( 



results 
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Proposition 1 

log (^|) = jf <&(*(«)) + ^(s)) 2 - g 2 (^)) 2 )cfe + jf g 2 (0(s))dW(s), 

(23) 

9(t) = 9(0)+ [ (q3(0(s))-q 2 (e(s))q 4 (d(s)))ds+ [ q 4 (9(s))dW(s), (24) 



(25) 



where 

qi(9) = an cos 2 9 + (a±2 + a 2 i) cos 9 sin 9 + a 22 sin 2 9, 
<? 2 (#) = fen cos 2 9 + (fei2 + fe 2 i) cos 9 sin + fe 22 sin 2 6>, 
q , 3 (6') = a 2 i cos 2 6* + (a 22 — 011) cos sin 9 — ai 2 sin 2 0, 
qi(9) = fe 2 i cos 2 + (fe 22 — fen) cos 9 sin 9 — fe i2 sin 2 9. 

As the expectation of the ltd stochastic integral is null, 

E [ t q 2 (9(s))dW(s) = 0, 
Jo 

the Lyapunov exponent is given by 

X = lim - t log (^|) = lim -E j\ qi {9{s)) + \{qMs)) 2 ~ 
Applying the Oseledec theorem, ifr(t) is ergodic, results that 

A = j [qM + ^(g 4 (0) 2 - q 2 (9))]p(9)d9, (26) 

where p(9) is the probability distribution of the process 9. □ 

An approximation of this distribution is calculated by solving the Fokker- 
Planck equation. Associated with equation f[2"4l for p = p(t, 9) we get 

% + lf9 {qm " q ^ q ^ ~ \w {q ' {e)2p) = °- (27) 

From (1271) results that the solution p(6 l ) of the Fokker-Planck equation is 
the solution of the following first order equation 

(-<&(*) + ?i W?4(f) + <&(% 5 (0))P(0) + i^WVW = Po, (28) 
where = -J and 

<?s(#) = -(612 + fe 2 i)sin20- (fe 22 - fen) cos 20. (29) 
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Proposition 2 If q±(9) 7^ 0, the solution of equation (|28p is given by 

pW = » (1+ 'f fl( " )fc) ' (30) 

where K is determined by the normality condition 



2- 



TTie function D is given by 



p(6)d6 = 1, (31) 
D(2n) - 1 



f* n D(u)du 



(32) 



□ 

A numerical solution of the phase distribution could be performed by a 
simple backward difference scheme. 
Let iV6l + and h = § . Let 

qi(i) = au cos 2 (ih) + (a i2 + a 2 i) cos(ih) sm(ih) + a 22 sin 2 (i/i), 

g 2 (i) = 6n cos 2 (ih) + (bi2 + b 2 i) cos(ih) sin(ih) + 6 22 sin 2 (ih), 

q 3 (i) = a 2 i cos 2 (ih) + (a 22 — an) cos(ih) sin(ih) — ai 2 sin 2 (i/i), (34) 

04(2) = b 2 i cos 2 (ih) + (b 22 — bu) cos(ih) sm(ih) — bi 2 sm 2 (ih), 

q 5 (i) = -(6 12 + b 2 i) sin(2ih) - (6 22 - b n ) cos(2ih). 

The p(i), i = 0, N is given by the following relations 

m = wo) + ^-'V w, 



where 



2/i 



2h(-q 3 (i) + q 2 (i)qS) + qi{f)qs{i)) + qS) 2 
The Lyapunov function is A = X(N), where 

N 



KN) = £>(i) + \{qSf - q2(i) 2 ))p(i)h. 



i=i 
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Proposition 3 // the matrix B is given by 

&n = a, bi2 = ~P, 621 = P, b 2 2 = a, 
probability distribution p(6) is given by 

p ^ = ~jj2 ex P{-^2~(( a ' 21 ~ fll2 ~~ + \( ai1 ~ a22 > cos2 ® + \( a21 ~~ 0,121 sin261 }' 



K 



P 2 



jo* ex P{^2 ((«2i - 012 - ocP)6 + |(an _ a 22) COS26 1 + |(a 2 i - ai 2 ) sin 20} 
and i/ie Lyapunov exponent is given by 



-d9, 



1 11 

A = -(an + «22 + P 2 - a 2 ) + -(an - ^22)02 + -(a 2 i + ai 2 )s 2 , 



where 



2- 



2- 



□ 



c 2 = J cos(26)p(9)d0, s 2 = J sm(2B)p(0)d9. 


A2 The Euler scheme. 

In general 2- dimensional case, the Euler scheme has the form: 



Xi(n + 1) = Xi(n) + fi(x(n))h + gi(x(n))Gi(n), i = 1, 2, (35) 
with Wiener process increment 

Gi(n) = Wi{{n + l)h) - W^nh), n = 0, N - 1, i = 1, 2, 

and £(n) = £(nh,uj), Gi(n) are generated using boxmuller method. 

It is shown that Euler scheme has the order for weak convergence 1, for 
sufficiently regular drift and diffusion coefficients. 

We assume that fi and a; in (131)]) are sufficiently smooth such that the 
following schemas are well defined. 

The second order Euler scheme is defined by the relations 

Q GAn) 2 - h 

Xi(n + l) = Xi(n) + fi(x{n))h + gi(x(n))Gi(n) + g i (x(n))—-—g i (x{n))— — 



+ 



fi(x(n)) 



dfi(x(n)) 
dxi (n) 



dxi{n) 

1/ / / \ \ 2 d 2 fi(x(n)) ih 2 

+ 2 (sw n )) ~ > <~ — — + 



dxi(n)dxi(n) 



9i{x{n)) 



df t {x{n)) 
dxi (n) 



%(x(n)) 1 2 d 2 gj(x(n)) -]hGj{n) 

+ fi(x(n)) +-{g i {x{n)) 

dxi(n) 2 0Xi(n)0Xi(n). 



,i = l,2, 
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where we used the random variables G^(n), i — 1,2. In [TT], it is shown that 
these schemes converge weakly with order 2. 
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